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I. INTRODUCTION 

The lattice Boltzmann (LB) method is a powerful ap- 
proach to hydrodynamics, with applications ranging from 
large Reynolds number flows to flows at a micron scale, 
porous media and multi-phase flows pj . The LB method 
solves a fully discrete kinetic equation for populations 
fi(x,t), designed in a way that it reproduces the tar- 
get equations of continuum mechanics in the hydrody- 
namic limit. Populations correspond to discrete veloc- 
ities Ci, i = 1,...,N, which fit into a regular spatial 
lattice with the nodes x. This enables a simple and ef- 
ficient 'stream-along-links-and-equilibrate-at-nodes' real- 
ization of the LB algorithm. 

In the case of incompressible flows, where the target 
equations are Navier-Stokes equations at low Mach num- 
ber, the LB method proves to be competitive to conven- 
tional methods of computational fluid dynamics [2|. In 
that case, the LB models on the so-called standard lat- 
tices with a relatively small number of velocities (N = 9 
in two dimensions, see Fig.[IJ and N = 15, 19, 27 in three) 
are available and most commonly used. In this paper, we 
will follow the usual nomenclature and indicate the mod- 
els as DMQN where M = 2, 3 is the spatial dimension 
and N the number of the discrete velocities of the model. 

However, situation is quite different with LB models 
for compressible thermal flows. In spite of a number of 
recent suggestions, a commonly accepted LB model for 
thermal flows has not yet been established, to the best 
of our knowledge. A reason for such a difference between 
the isothermal and thermal cases is mainly because it is 
not straightforward to incorporate the temperature into 
the lattice equilibrium when using the standard lattices, 
and simultaneously to satisfy a number of conditions for 
recovering the Navier-Stokes and Fourier equations for 
compressible flows Q. At present, two major approaches 
to constructing thermal LB models can be distinguished. 
In the first approach, lattices with a larger number of dis- 
crete velocities or off-lattice velocity sets are considered 
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to enable local energy conservation and isotropy [J, |5j . In 
the second approach, two copies of the standard lattices 
are considered, one of which is designed to treat the den- 
sity and momentum and the other - the energy density. 
The standard isothermal LB model on the first lattice 
is considered, and the coupling between the lattices is 
enhanced by introducing force terms in order to recover 
viscous heat dissipation However, both these ap- 
proaches inevitably deal with more discrete speeds than 
the standard single-lattice LB models which leads to re- 
ducing efficiency. 

Recently, a so-called consistent LB model has been in- 
troduced in Ref. 0- It has been shown that it is possi- 
ble to construct the LB model with energy conservation 
on the standard lattices in such a way that the spurious 
bulk viscosity of isothermal LB models is eliminated, and 
that the Navier-Stokes and Fourier equations are recov- 
ered once the temperature is varied in a small neighbor- 
hood of a reference temperature. However, the consis- 
tent LB model is limited to weakly compressible flows, 
and is not yet sufficient to simulate thermal flows. The 
reason for that is the low symmetry of the standard lat- 
tices which leads to accumulation of deviations of the re- 
sulting hydrodynamic equations from their Navier-Stokes 
and Fourier counterparts. It should be noted, however, 
that the consistent LB model with energy conservation 
on standard lattices should be considered as a better and 
more natural starting point for a development of the ther- 
mal LB model as compared to the isothermal model. 

In this paper, we complete the program of construct- 
ing the thermal LB model on standard lattices as an ex- 
tension of the consistent LB model A new thermal 
lattice Boltzmann model is derived on the most com- 
monly used D2Q9 lattice as follows: In section UH wc 
revisit the derivation of the local equilibrium of the ther- 
mal model, using the method of guided equilibrium [f|. 
This enables to enhance Galilean invariance of the con- 
sistent LB method for large temperature variations. We 
shall also identify the remaining deviations in the higher- 
order moments due to the lattice constraints. In section 
IIII1 the impact of these deviations on the hydrodynamic 
equations is identified (algebraic details of the deriva- 
tion are presented in the Appendix. In section UVl addi- 
tional terms are introduced in the Boltzmann equation in 
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FIG. 1: The D2Q9 velocity set. 



such a way that all the deviations of the hydrodynamic 
equations from their Navier-Stokes and Fourier counter- 
parts are eliminated. In addition, further terms are in- 
troduced to represent external forces, as well as to tune 
the Prandtl number of the model. In section [V] the dis- 
cretization in time and space of the kinetic equation of 
section IIVI is performed, and a simple numerical algo- 
rithm is outlined. In section IVH numerical examples are 
provided. Simulations of the thermal Couette flow and 
of the Rayleigh-Benard natural convection show excel- 
lent agreement with theoretical results. The numerical 
implementation is summarized in a compact form in sec- 
tion IVIII Finally, results are discussed in section IVIII1 



II. LOCAL EQUILIBRIUM 



A. Consistent lattice Boltzmann method 



We begin with reminding the construction of the local 
equilibrium in the consistent lattice Boltzmann method 
0, Q on the planar square lattice with nine velocities Ci a , 
i = 0,...,8 (see Fig. EE): 



c x = {0,1,0,-1,0,1,-1,-1,1} 
Cy = {0,0,1,0,-1,1,1,-1,-1} 



(1) 



In the consistent lattice Boltzmann method Q , the local 
conservations include mass, momentum and energy, 
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where p, j a — pu a and T are the density, momentum and 
temperature fields, correspondingly. Note that the con- 
sistent lattice Boltzmann construction includes energy 
conservation among the constraints (H]). Note that this 
is at variance with the standard isothermal lattice Boltz- 
mann method on the same lattice [ll[ where the entropy 
function is minimized under the constraints of mass and 
momentum only [lfj. 

Derivation of the equilibrium populations as a mini- 
mizer of the entropy function ^ under constraints ^ 
can be found in Refs. @, Q, and is not reproduced 
here. We remind that at the reference temperature 
To = 1/3, the consistent lattice Boltzmann model recov- 
ers the Navier-Stokes-Fourier hydrodynamic equations 
and results in the nonequilibrium pressure tensor without 
a spurious bulk viscosity of the standard lattice Boltz- 
mann model. However, the consistent lattice Boltzmann 
model itself is not yet capable of simulating accurately 
thermal flows with significant temperature and density 
variations. In particular, the equilibrium pressure ten- 
sor of the consistent lattice Boltzmann model shows de- 
viations of the order of u 2 AT (where AT is a devi- 
ation from the reference temperature) from the target 
Maxwcll-Boltzmann form, 



P^ = pT5 a , 



JaJ/3 



(5) 



Hence, as the first step towards establishing the lat- 
tice Boltzmann method for thermal flow simulations, we 
need to modify the construction of the equilibrium which 
would remove this deviation from the pressure tensor. 



Equilibrium populations are obtained from minimization 
of the entropy function H under constraints provided by 
local conservation laws. For the D2Q9 model, the en- 
tropy function has the form [1 01 ] . 

i=0 v 7 

where the weights Wi are: 

W= ^ {16, 4, 4, 4, 4, 1,1, 1,1}. (3) 



B. Guided local equilibrium 

In order to remove the aforementioned deviation in 
the equilibrium pressure tensor, we use the method of 
guided equilibrium introduced in Q for a generic lattice 
Boltzmann model. Following Q, we minimize the en- 
tropy function ^ under an extended set of conditions 
which includes the local conservations 0$ and the condi- 
tion which stipulates that the equilibrium pressure tensor 
is in the MB form ©, 

8 . . 

^ CiaCipfP = pT8 a p + --V-.. (6) 
i=0 " 



3 



Minimization of the entropy function ([2]) under the ex- 
tended list of constraints, ([U and © leads to the equi- 
librium populations of the form 



n 



j exp{/J + C, a C ia + Xaf3C ia Ci0 + 7C?}, (7) 



where summation convention in two repeated spatial in- 
dices is assumed, and where p, £ Q , \ap, 7 are Lagrange 
multipliers of the constraints. Their values are found 
upon substitution of into (0J and ([5]). Equilibrium 
at zero velocity (j a = 0) can be found exactly and coin- 
cides with the one reported in 0, Q . Equilibria at non- 
zero velocity are then derived by perturbation around the 
zero-velocity state, and /? q are represented in terms of a 
series in the momentum. For what will follow, expansion 
up to the fourth order in velocity is required. This so- 
lution procedure is quite standard, hence, we reproduce 
only the result in terms of the velocity u a = j a /p: 
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Note that equilibrium populations at zero velocity re- 
main positive for temperature values < T < 1. The 
equilibrium pressure tensor satisfies the MB relation per 
construction, 



poq _ p MB 
r u8 — r aB ■ 



(9) 



The properties of the higher-order moments in the equi- 
librium HO) are studied in the next section. 



Using the equilibrium populations ([8]) to evaluate the 
functions 



Qap-y ~ Ci a Cif3Ci 7 f i q , 
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and introducing, for a generic moment M, a devia- 
tion M' of its equilibrium value M eq from the Maxwell- 
Boltzmann value Af MB , 



M eq = m mb , M / 



(12) 



we find the following deviations of the third- and fourth- 
order moments: 
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Note that the off-diagonal elements of the third-order mo- 
ment Q^jL, already satisfy the Maxwell-Boltzmann rela- 
tions, 



jeq _ amb n e q — n MB (ia\ 

z XX y *>ixxy> ^Cyy X ^Cyy X ' \ J 



Thus, the deviation of the contracted third-order mo- 
ment, ql q = £)i=o c ia cf ft q (the energy flux) is due to 
the deviation of the diagonal elements, 

q' a = Q' aaa - (is) 



C. Deviations in higher-order moments 

Apart from the equilibrium pressure tensor, also the 
third- and fourth-order moments of the equilibrium must 
satisfy the Maxwell-Boltzmann relations in order to re- 
cover the Navier-Stokes and the temperature equation in 
the hydrodynamic limit, 
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Deviations of the diagonal components of the third- 
order moments are well known for the low-symmetry 
D2Q9 lattice, and stem from the fact that the velocity 
set satisfies a relation, cf Q = c$ a . This lattice constraint 
precludes construction of equilibrium different from (JSJ 
in such a way that also the energy flux would be guided 
by the Maxwell-Boltzmann relation. 

In the next section, we shall identify the impact of the 
deviations (|13|) on the hydrodynamic equations. Once 
this will be done, we shall remove the anomalous terms 
in the hydrodynamic equations by introducing correction 
terms into the kinetic equation. 
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III. DEVIATIONS IN THE HYDRODYNAMIC 
EQUATIONS 



B. Deviations in the momentum equation 



A. Bear BGK model 

In this section, the impact of deviations of the moments 
(|13|) on the hydrodynamic equations will be identified 
using the single relaxation time Bhatnagar-Gross-Krook 
(BGK) model for the collision integral: 



dtfi + c la d a f t = —(fi- fP), 

T 



(16) 



To the first order of the Chapman-Enskog expansion, 
the momentum equation reads: 



d t j a + dpi*? + dpPlf = -d y P" , (17) 



where -P"^ q is the nonequilibrium pressure tensor in the 
form required for the Navier-Stokes equation, 



where r is the relaxation time, and the equilibrium pop- 
ulations are given by Eq. p]). The hydrodynamic limit 
of the BGK kinetic equation (JTSJ) can be studied via the 
Chapman-Enskog expansion. Details of the Chapman- 
Enskog analysis are given in the Appendix [SJ Here we 
shall summarize the results of this analysis. 
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while P'^a is the deviation of the nonequilibrium pressure 
tensor from P"^ q : 
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In order to recover the Navier-Stokes equation, we will 
need to introduce additional terms into Eq. (fi7j|) in or- 
der to annihilate the right hand side of the momentum 
equation (TT7|) (see section ITV)) . 



C. Deviations in the energy equation 



Similarly, the energy density equation reads, 
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is the nonequilibrium energy flux required in the Fourier 
equation for the energy density. The first term in (|21D is 
the Fourier law of energy dissipation, whereas the second 
term represents viscous dissipation. Again, terms in the 
right hand side of (|20|) express the deviation of the energy 
equation from the required form. The first of these terms 
is given by Eqs. (fTS)) and (fT3")) , while the result for the 
second term, q^, is given by the following expression: 
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It should be noted that, in certain situations, many terms 
in (j2"2")l and (|2"3"|) can be safely neglected. For example, 
deviations of the order j 4 affect the viscous heat dissipa- 
tion and can be ignored in cases where viscous heating 



5 



is unimportant. However, we shall use the full expres- 
sions (|2"2")) and ([25]) for the deviation in the numerical 
implementation below. 

The deviations in the hydrodynamic equations identi- 
fied in this section are due to the lattice constraints of 
the D2Q9 lattice, and they cannot be removed by in- 
troducing a set of the equilibrium distribution functions 
different from ijHJ) . These deviations can be removed only 
by introducing correction terms into the kinetic equation 
which we address in the next section. 



IV. CORRECTION OF BGK EQUATION 

To this end, we have computed the deviations due to 
the lattice constraints. In this section, an efficient way 
to remove these deviations by adding correction terms in 
the Boltzmann equation is introduced: 



Among these equations only nine are independent, and 
the unique solution is: 
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dtfi + c ia d a fi =—(fi- O + 9 t + *i. (24) 
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The purpose of the terms is to correct the momentum 
equation, while the purpose of the $^ terms is to correct 
the energy equation. 



A. Momentum equation correction terms 9i 

We require that the 'J'i term affects only the momen- 
tum equation and delivers there a term which compen- 
sates — d 1 P 1 '. Thus, 



B. Energy equation correction terms $i 



Similarly, it is required that <&i terms appear only in 
the energy equation, 
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The solution of that system is: 
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C. Variable Prandtl number 

A thermal lattice Boltzman model should be able to 
simulate fluids with an arbitrary Prandtl number. The 
Prandtl number Pr is defined as the ratio of the viscosity 
and the thermal conductivity, 

Pr=^, (29) 

K 

where c p is specific heat under a constant pressure. We 
here suggest the following simple way to change the 
Prandtl number in the present model by adjusting the 
correction term In order to do this, we first note 
that, with the use of the guided equilibrium populations, 
the BGK model (TT^|) leads to Pr = 4, same as the orig- 
inal consistent LB method [9(. Applying the correction 
for the energy equation mentioned in the previous sec- 
tion, the Prandtl number becomes Pr = 1, as in the 



standard continuous BGK model. This happens because 
the correction compensates, among the others, the term, 

3rpTd a T, 

(first term in (I22]) ). responsible for the change of the ther- 
mal conductivity of the model. Thus, arbitrary values of 
Pr can be obtained by applying the following transfor- 
mation in (|22p and thereby altering the counter-term $^ 
in is term: 

(4 - Pr) 

3rpTd a T -> ^ -rpTd a T. (30) 

Pr 

Thus, with the simple transform (1301) in the energy cor- 
rection term kinetic equation (j2"4")) recovers the hydro- 
dynamic equations with a predetermined Prandtl num- 
ber. In section IVI( we shall validate this in the thermal 
Couette flow for fluids with various Prandtl numbers. 
D. Gravity force 

Flows subject to external forces, such as gravity, are a 
major concern in many applications. In order to incor- 
porate the force pg, where g is the acceleration due to 
gravity, it suffices to apply the following transformation 
when calculating the correction terms tyf. 

dpKp - dpP*p + P9 a - (31) 

Any additional physics can be incorporated into the 
present model by similar considerations. 



E. Corrected hydrodynamic equations 

Taking into account all the aforementioned correc- 
tion terms and <&i in the kinetic equation (|24|) . and 
using again the Chapman-Enskog analysis, we obtain 
the two-dimensional Navier-Stokes and Fourier hydrody- 
namic equations for density, momentum and tempera- 
ture: 



dtp = 
d t u a = 
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(32) 



The fluid corresponds to the ideal gas equation of state, 
p = pT, with specific heats c v = 1.0 and c p — 2.0, in lat- 
tice units, and with the adiabatic exponent 7 = c p /c v = 
2.0. The viscosity coefficient p and thermal conductivity 



k can be identified as 



p = rpT, k = ^-rpT. 

Pr 



(33) 
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The Prandtl number is a tunable parameter and can take 
arbitrary values. A wide range of gases and fluids whose 
physics is described by equations (|32|) can be simulated. 
We now proceed with a time-space discretization of the 
kinetic equation (|24]) . 



Ref. [6j. First, we integrate over the time step 5t along 
the characteristics: 



V. LATTICE BOLTZMANN DISCRETIZATION 



Derivation of the lattice Boltzmann scheme for the ki- 
netic equation (f24f proceeds essentially along the lines of 



t+st 
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fi(x + CiSt,t + 6t) = fi(x,t) + 
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W))dt' 



(34) 



The time integrals of the collision term as well as of the 
correction terms is evaluated using the second-order ac- 
curate trapezoidal rule. Second, in order to establish a 
semi-implicit scheme, the following transformation is ap- 
plied @: 



fi + 2 T 



/D -?[*< + *<] (35) 



This leads to a semi-explicit scheme, 



9t+st = 9t + 
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(36) 

Note that the simple scheme (|36|) is semi-explicit due to 
the presence of the correction terms \&j and <&i (and not 
fully explicit, as in the standard case without any correc- 
tions). The scheme utilizes the transformed populations 
g. The equilibrium / oq can be computed using equations 
that relate the locally conserved moments of the popu- 
lations / with those of the ^-populations. In order to 
do this, we evaluate the moments of the transformation 
(El: 
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jyU) = jy(9) + - Ci v 
T(f) = 
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The above set of equations can be simplified using Eqs. 
(HI and HMD: 



P(f) 
Mf) 

Jvif) 
T(f) = 



p(g), 



st 



■2 




Finally, discretization in space is done as in the standard 
lattice Boltzmann method: if x is a grid node then also 
x ± CiSt is the grid node. 

In the simulation, the following algorithm was imple- 
mented for the collision step: 

Step 1. Calculate p, j a , T using gUS]), (^P^)* -1 , 
(d a (q' a + q'a))*' 1 and g t values; 

Step 2. Calculate (9 7 P^ 7 )*, (d a (q' a + g"))* using values of 
p, j a , T from Step 1; 




j 2 (f) , st 



?E^ t (/) V40) 



Step 3. Calculate again p, j a , T using (|4HI44[) . g t and the 
values calculated in Step 2; 

Step 4. Use p, j a , T from Step 3 for the calculation of the 
equilibrium values (jHJ); 

Step 5. Use (SyP^)*, {d a {q' a + q'^Y from Step 2 in the 
discrete equation (|36[) along with the equilibrium 
values calculated in Step 4. 

The terms 9 7 P^' 7 and d a {q' a + q'£) are evaluated using a 
second-order accurate finite-difference scheme. 

A few comments on the computational efficiency of our 
model are in order. First, the present thermal model is 
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only about 2.7 times slower as compared to the stan- 
dard isothermal LB method on the same D2Q9 lattice 
due to evaluations of the deviation terms (fT9|) and (f22|) . 
Note that, depending on the desired accuracy and the 
optimization level, faster and simpler algorithms can be 
readily designed. For example, for the natural convec- 
tion problem considered below in section IVI1 deviations 
of the order j 3 , j 2 AT, and j 4 can be safely neglected in 
the energy equation, which removes the computational 
overhead with respect to the standard isothermal D2Q9 
code without sacrificing the accuracy of the results in the 
thermal case. It should be stressed that optimization of 
the code was not pursuit in the present study so that the 
above estimates are only qualitative. 

Second, as compared to the thermal model on two lat- 
tices @ in terms of memory use, the present model re- 
quires storage of twelve values per node, that is, nine for 
the populations gt, two for the terms (9 7 P" ) , and one 

for the term {d a (q' a +<;"))*■ In the case of the two dis- 
tribution function thermal models @, the storage of at 
least eighteen values per node is required, without taking 
into account the additional computational effort needed 
in order to recover viscous dissipation terms. 

These estimates show that the present one-lattice ther- 
mal model is viable. We shall now proceed with a nu- 
merical validation of the present scheme. 



VI. NUMERICAL EXAMPLES 

In this section, the thermal model developed above 
is validated numerically. Equilibrium populations are 
expanded till the fourth order in the velocity and the 
corrections described in previous sections are applied. 
The same model is used in all simulations, even though 
viscous heat dissipation could have been ignored in the 
Rayleigh-Benard setup. 



where P = P cq = I 
tensor and N = P TT — P, 



yy 



^yy is the trace of the pressure 
is the normal stress difference. 



On one hand, using the guided equilibrium © in (|4"5|) 
results in the constant pressure in the whole domain: 



pT = const. 



(46) 



This is consistent with the correct Maxwell-Boltzmann 
relation for the pressure tensor. On the other hand, as we 
have already mentioned in section[TTl the equilibrium nor- 
mal stress difference N cq of the consistent LB model @, Q 
exhibits a deviation of the order ~ j 2 AT, where AT is a 
variation around the reference temperature To = 1/3. In 
this case, equation (|45|) yields a different result, namely, 
the pressure p = pT varies along the y-axis according to 
the variation of the momentum, 



P T-^1^1 = const. 



(47) 
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In Fig. [21 simulation results for the pressure in the Cou 
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A. Necessity of guided equilibrium 

In our first example, we compare the three different 
models: the original consistent LB model 7, 9], the 
model using the guided equilibrium ((SJ) without any cor- 
rection terms, and the full model after the correction 
terms are applied. In order to demonstrate the ne- 
cessity of the guided equilibrium, we simulate the two- 
dimensional Couette flow between two parallel moving 
plates at different temperatures. Isothermal walls are 
separated by a distance H and move parallel to the x- 
axis. The diffusive wall boundary conditions [l| is imple- 
mented for the isothermal walls, and periodic boundary 
condition is applied for the rest. Analysis similar to the 
one presented in [l3| results in the following relation at 
the steady state: 

P — N cq = const, (45) 



FIG. 2: Pressure p = pT as a function of the reduced y coor- 
dinate (ynorm = y/H, H is the distance between the plates) 
in the thermal Couette flow. Line: analytical solution (146 [) : 
Open circles: LB model with guided equilibrium; Filled cir- 
cles: LB model with guided equilibrium and correction terms. 
Diamonds: LB model 

UM- Reynolds number Re = 320. 

ette flow with all the three models are presented and 
compared with the analytical solution (f46|) . While the 
variation of the pressure p — pT away from the constant 
value is numerically small for the original model 0, Q , it 
is still visible, and it verifies the analytical solution ([47]) 
rather than (|46|) . On the other hand, the LB model with 
the guided equilibrium verifies the constancy of the pres- 
sure with machine precision. Same holds also for the LB 
model with the guided equilibrium when all the correc- 
tion terms are applied, as expected. This demonstrates 
necessity of the guided equilibrium in our construction. 
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B. Viscous heat dissipation 



We simulate the thermal Couette flow as in the pre- 
ceding section but for a small temperature differences. In 
this case, the viscous heat dissipation is important and 
affects the temperature profile. As is well known, the 
Navier-Stokes and Fourier equations (|32| predict that the 
temperature profile depends on the Prandtl number, Pr, 
and on the Eckert number, Ec, where 



Ec = 



CpAT' 



(48) 



with u the difference of the velocities of the plates, and 
AT the difference of temperatures of the plates. This 
simulation enables to verify that the correction terms do 
recover the right Prandtl number. In the simulation, we 
used Ec = 3 and AT = 2xlO" 4 T , where T = 1/3 is 
the reference temperature. Three different cases were 
considered: 

Case 1. LB with the guided equilibrium ([5]) and the cor- 
rections Vfi and $i which deliver Pr = 1 (that is, 
without the transformation (130|) : 



Case 2. As in Case 1 plus the transformation (|30|) to match 
Pr = 0.71 (air). 

Case 2. As in Case 1 plus the transformation ([3171) to match 
Pr = 4. 

Simulation results for the temperature profile are pre- 
sented in Fig. [3J and are in excellent agreement with 
the analytical solution. For the current setup, with the 
present implementation of boundary condition, simula- 
tion results remain in agreement with the analytical so- 
lution for Prandtl numbers ranging from 0.01 to 20. 



C. Rayleigh-Benard convection 

The Rayleigh-Benard convection flow is a classical 
benchmark on the thermal models. The fluid is enclosed 
between two parallel stationary walls, the hot (bottom) 
and the cold (top), and experiences the gravity force. 
Density variations caused by the temperature variations 
drive the flow, while the viscosity will counteract to equi- 
librate it. 

In our LB model, gravity is implemented using the cor- 
rection (|3"Tj) . We operate the model in the weakly com- 
pressible regime, however, without using the Boussinesq 
approximation. For that we set a temperature difference 
of the order of 3% of the reference temperature To = 1/3. 
At the top and the bottom walls we apply the diffusive 
wall boundary condition [l2j , whereas periodic boundary 
condition is applied on the vertical walls. For the nodes 
that belong to the isothermal walls, gravitational force 
was not applied. 




0.4 0.6 

y-direction 



FIG. 3: Temperature profile in the thermal Couette 
flow. Plotted is the reduced temperature, Xkorm — (T — 
T co id)/(Thot — T co id), versus the reduced y-coordinate (as in 
Fig-El- Eckert number Ec = 3. Line: Analytical solution for 
Pr = 4, Pr = 1, and Pr = 0.71. Reynolds number Re = 200. 
Symbol: Simulation results for the three cases (see text). Di- 
amonds: Case 1; Squares: Case 2 (air); Circles: Case 3. 



The Prandtl number used corresponds to air, Pr = 
0.71. The dimensionless number that characterizes the 
flow is the Rayleigh number Ra, 



Ra 



Pr gPATH 3 



(49) 



where g is the gravity acceleration, AT is the wall tem- 
perature difference, H is the distance between the walls 
and v is the kinematic viscosity of the fluid. For the ideal 
gas equation of state, the thermal expansion coefficient 
(3 is defined as: 



/* = -- i£ 



I ( dp_ 

P \dT 



1 

T' 



(50) 



For the computation of the Rayleigh number, we used 
the thermal expansion coefficient evaluated at the refer- 
ence temperature T) = 1/3, that is, (3 — 3. The heat 
transfer is described by the Nusselt number, Nu, defined 
as the ratio between convective heat transport to the heat 
transport due to temperature conduction: 



Nu = 1 



(UyT)H 

kAT 



(51) 



Here (u y T) denotes the average over the convection layer 
and k is the thermal conductivity of the model (|33[) . For 
the current simulations, a computation domain with the 
aspect ratio 2 : 1 was considered. 

For this setup, the critical Rayleigh number was found 
to be Ra cr = 1700 ± 10, which is consistent with the 
reference data for the Boussinesq approximation [3, [H| • 
Any flow field was dissipated for values less than Ra cr . 
On the contrary, for values larger than Ra cr , the flow 
develops two or more cells (vortices) depending on the 
initial conditions. 
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FIG. 4: Contour plot with 19 equidistant iso-temperature 
lines for Ra = 10 4 (top) and Ra = 5xl0 4 (bottom). 
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FIG. 5: Stream function contours for Ra = 10 4 (top) and 
Ra = 5xl0 4 (bottom). Difference between contour lines is 

0.025 in Units Of pmaxUmaxH. 

For different Rayleigh numbers, simulations for the 
51, 101 and 201 grid nodes in the y-direction were per- 
formed. Extrapolating the obtained values of Nusselt 
number, an estimate of the final converged solution can 
be done (Nur). In Fig. [4j the isotherms of Ra = 10 4 and 
Ra = 5x10 are plotted for the case of 101 grid nodes 
in the y-direction. In Fig. [5l the contours of the stream 



function of the compressible flow field for Ra = 10 4 and 
Ra = 5x1 4 are plotted. The extrapolated converged val- 
ues of Nusselt number at various Rayleigh numbers are 




i *-* , , , , , , 

10 3 10 4 

Rayleigh number 



FIG. 6: Rayleigh Benard convective flow. Nusselt number 
vs. Rayleigh number. Diamonds: The current LB model; 
Triangles: Reference data of Ref. [l4|; Line: Empirical power- 
law Nu = 1.56(Ra/Ra c ) - 296 [3. 




y-direction grid nodes 



FIG. 7: Grid convergence study for Ra = 3xl0 4 . Symbol: 
The relative error between the simulation result Nu and the 
converged solution Nur,. Line is a fitted curve which reveals 
second order convergence. 



plotted in Fig [61 and compared with the standard refer- 
ence data pj], LL5| in Table Q] The present model is found 
to be in good agreement with [3| , and it also agrees well 
with the thermal LB models on double lattices where 
the temperature dynamics is treated as a passive scalar 
[UEB]. It should be noted that for Ra up to 10 4 a uniform 
grid with 51 nodes in the y-direction was sufficient, while 
for larger values of Ra, larger grids provide more accurate 
results. In Fig. [Jj a grid convergence study is performed. 
Natural convection of Ra = 3xl0 4 in the same setup is 
considered. Fig. [7j reveals the second-order accuracy of 
the numerical scheme. 
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Rayleigh number 


Nusselt number (Nur) 


Difference from Ref. [141 


2500 


1.474 


0.07 % 


5000 


2.104 


0.57 % 


10000 


2.644 


0.64 % 


30000 


3.605 


1.56 % 


50000 


4.133 


2.64 % 



TABLE I: Simulation results of the present model compared 
to the simulation results of Ref. [l4j| . Left column: Rayleigh 
number; Central column: Nusselt number; Right column: Dif- 
ference in percent from the values of Ref. [1J| . 




x-direction 

FIG. 8: (Color online) Snapshot of the temperature field at 
Ra = 10 8 , Pr = 0.71, using a uniform grid of 51 nodes for 
the y-direction. Contours of reduced temperature, T norm = 
(T - T cold )/(T hot - T cold ) are plotted. 

Finally, it should be mentioned that on a grid with only 
51 nodes in the y-direction, the code was run stably for 
Ra as high as at least 10 8 which proves a good numerical 
stability of the algorithm. In Fig. [5J a snapshot of the 
temperature field is presented at Ra = 10 8 . However, a 
study of the natural convection at high Rayleigh numbers 
is out of the scope of this paper. 



VII. SUMMARY 

A reader who wishes to implement this model should 
execute the following steps: 

• Use guided equilibrium populations (|S]); 

• Calculate the correction terms <Fi and with a 
second-order finite-difference scheme. 

• Use the discretized equation (136)) . 

VIII. CONCLUSION 

In this paper, we have developed a new lattice Boltz- 
mann model for simulation of thermal flows with the 



standard and most commonly used D2Q9 lattice. The 
important starting point of the derivation is the consis- 
tent LB method p}. Unlike the previous approaches, the 
consistent LB model on standard lattices already includes 
energy as a locally conserved quantity, the nonequilib- 
rium stress tensor is free of a spurious bulk viscosity, 
and a part of the moment relations satisfies the required 
Maxwell-Boltzmann relations. Hence, by using the con- 
sistent LB model, it becomes unnecessary to introduce a 
separate lattice for the energy field. 

We have considerably refined the consistent LB model 
in two steps. First, following the concept of guided equi- 
librium [8| , we recovered Galilean invariance of the pres- 
sure tensor at non-vanishing variations of the tempera- 
ture. Second, we have identified the remaining devia- 
tions in the Navier-Stokes-Fourier equations, and have 
removed them by adding compensating terms into the 
kinetic equation. This step allowed us also to make the 
Prandtl number a tunable parameter, and to add exter- 
nal forces. It should be noted that additional terms also 
appear in the method based on two lattices [|| but there 
is no relation between the both. On the other hand, a 
numerical implementation of the kinetic equation with 
additional terms here has practically not much differ- 
ent from the method first developed in [||. Finally, it 
is straightforward to extend the present model to three 
dimensions based on the D3Q27 consistent LB model 0] ■ 

We shall conclude this paper with a general comment 
on the construction of thermal lattice Boltzmann models. 
As we already mentioned in the introduction, there are 
two basic directions for such a development, one which 
uses larger and more isotropic sets of velocities, and the 
other which uses a small number of velocities and intro- 
duces non-local (dependent on the gradients of the fields) 
corrections to compensate for the deficit of lattice sym- 
metry. This situation resembles the classical dilemma 
of extending the hydrodynamics into beyond the Navier- 
Stokes level: One way is the higher-order hydrodynam- 
ics (Burnett or super-Burnett) which means higher-order 
derivatives and more non-locality in the equations for the 
standard hydrodynamic fields. The other is to take more 
moments while staying local (Grad's moment method). 
Both approaches have positive and negative aspects. In 
our case, the attractive features of large velocity sets 
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(locality, isotropy, possibility of the exact treatment of 
propagation) should be waged against the necessity of 
handling many more fields (populations). On the other 
hand, the obvious attractive side of the small lattices ap- 
proach has to be opposed by the fact that one has to 
evaluate more terms dependent on the gradients of the 
fields. While this proves to be possible, still, one needs 
to be careful about numerical errors brought in by such 
procedures. We believe that both these approaches, at 
present, show there advantages and disadvantages, so the 
decision about which to prefer should be the focus of fu- 
ture studies. 

We are indebted to our colleagues S. Ansumali, S. Ar- 
cidiacono, K. Boulouchos, S. Chikatamarla, C. Frouza- 
kis and J. Mantzaras for their help, questions and dis- 
cussions. We gratefully acknowledge support by BFE 
Project No. 100862, CCEM-CH and by ETH Project 
No. 0-20235-05. 



On the first order in r we need equations for the locally 
conserved fields only, 



r j 2 



d\ >p 
u t Jot 

12 P T 



0, 



-d o (1) 



(A10) 
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Function is represented as a sum of two terms, -P^J 1 , 
which results in the correct Navier-Stokes term in the mo- 
mentum equation, and P^g, representing the deviation: 
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APPENDIX A: DERIVATION OF 
HYDRODYNAMIC EQUATIONS 



In order to derive Pf£j from (|A7[) . we note that the left 
hand side can be computed by chain rule, 



We apply the Chapman-Enskog expansion in order 
to derive hydrodynamic equations corresponding to the 
BGK equation (JTHJ) with the equilibrium §8§ Populations 
fi and the time derivative operator are expanded into 
powers of the relaxation parameter r 



fi 
di 



^n=0 



(Al) 
(A2) 



The non-equilibrium pressure tensor and energy flux to 
first order in r are defined as: 



Using 
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6j) and (|A16|) in (TATf . we obtain 

* (7)^(7) iJ 

(A18) 
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On the zeroth order, equations for density, momentum, 
pressure tensor and energy flux are: 
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(A6) 

(A7) 
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where Q a q g^ and R A are the third- and the fourth-order 
moments evaluated at the local equilibrium ([8]) (explicit 
form of these function was given in section [TTJ) . From 
(|A6|) . we derive the zeroth-order equation for the tem- 
perature, 



d W T = _ TQa J J* j _ J^ daT _ ±.d a C (A9) 



3c 



1 



Substituting (|A17|) and (|A18|) into (|A11|) . and combining 
the latter with (|A5[) . and also combining (|A4|) with (jAlO|) . 
we derive the equations for the density and velocity u a = 
ja/p to first order, 



d t p 
d t u a 



-d a {pu a ), 
-u B dgu a - 

--d P" 



-8 a { P T)- -d B n a0 
p p 



(A19) 



(A20) 



where 



-T~pT [d 8 u a + d a u g - (9 7 u 7 )(5 Q/3 ] , (A21) 



is the nonequilibrium pressure tensor of a Newtonian 
fluid, and P£ (|A18|) is the deviation. Expanding Eq. 
(|A18|) . we arrive at the explicit form of the deviation, 



Eq. |[1J). 

Derivation of the equation for the energy (or, equiva- 
lently, for the temperature) is done in a similar way upon 
computing the first-order correction from equation 
The resulting equation for the temperature reads: 
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dtT = —u a d a T — Td a u a 

[d a {2 T pTd a T) ~ {df3U a )n a p} 

-Tp^-Tp^ (A22) 

The last two terms in this equation represent the devi- 



ation. Note that, unlike in the case of the momentum 
equation (|A20[) . the deviation also includes a zero-order 
term q' K , This happens because, with the choice of the 
guided equilibrium |5]), we have set the equilibrium pres- 
sure tensor in the required Maxwell-Boltzmann form but 
not the equilibrium energy flux (why the latter is impos- 
sible on the D2Q9 was explained in section [TTJ) . Finally, 
the first-order deviation (last term in (|A22[) ) has the form 
given by Eqs. (|22) and (|23|) . 



[1] S. Succi, The Lattice Boltzmann Equation for Fluid Dy- 
namics and Beyond (Oxford University Press, Oxford, 
2001). 

[2] H. Chen, S. Kandasamy, S. Orszag, R. Shock, S. Succi, 

and V. Yakhot, Science 301, 633 (2003). 
[3] S. Succi, I. V. Karlin, and H. Chen, Rev. Mod. Phys. 74, 

1203 (2002). 

[4] S. Ansumali, I. V. Karlin, and H. C. Ottinger, Phys. Rev. 

Lett. 94, 080602 (2005). 
[5] S. S. Chikatamarla and I. V. Karlin, Phys. Rev. Lett. 97, 

190601 (2006). 

[6] X. He, S. Chen, and G. D. Doolen, J. Comput. Phys. 

146, 282 (1998). 
[7] S. Ansumali and I. V. Karlin, Phys. Rev. Lett. 95, 260605 

(2005). 

[8] I. V. Karlin and S. Succi, Phys. Rev. E 56, R4053 (1998). 



[9] N. I. Prasianakis, S. S. Chikatamarla, I. V. Karlin, S. An- 
sumali, and K. B. Boulouchos, Math, and Comp. in Sim- 
ulation 72, 179 (2006). 

[10] I. V. Karlin, A. Ferrante, and H. C. Ottinger, Europhys. 
Lett. 47, 182 (1999). 

[11] Y. H. Qian, D. d'Humieres, and P. Lallemand, Europhys. 
Lett. 17, 479 (1992). 

[12] S. Ansumali and I. V. Karlin, Phys. Rev. E 66, 026311 
(2002). 

[13] S. Ansumali, I. V. Karlin, S. Arcidiacono, A. Abbas, and 
N. Prasianakis, Phys. Rev. Lett. 98, 124502 (2007). 

[14] R. M. Clever and F. H. Busse, J. Fluid Mech. 65, 625 
(1974). 

[15] F. H. Busse, Rep. Prog. Phys. 41, 1929 (1978). 
[16] X. Shan, Phys. Rev. E. 55, 2780 (1997). 



